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If the lighter top squark has no kinematically allowed two-body decays that conserve flavor, 
then it will live long enough to form hadronic bound states. The observation of the diphoton 
decays of stoponium could then provide a uniquely precise measurement of the top squark 
mass. In this paper, we calculate the cross section for the production of stoponium in a 
hadron collider at next-to-leading order (NLO) in QCD. We present numerical results for 
the cross section for production of stoponium at the LHC and study the dependence on beam 
energy, stoponium mass, and the renormalization and factorization scale. The cross-section 
is substantially increased by the NLO corrections, counteracting a corresponding decrease 
found earlier in the NLO diphoton branching ratio. 
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I. INTRODUCTION 

In low-energy supersymmetry, the usual collider signatures feature missing energy carried away 
by stable, neutral, weakly-interacting particles. Therefore, there are typically no kinematic mass 
peaks that could be used to precisely measure the individual superpartner masses. Other kinematic 
features can be used to determine mass differences with some precision, if the decay chains are 
favorable, but the overall mass scale of the sparticle spectrum is likely to be much more difficult 
to pin down accurately. 

One exception to this may occur if new particles can form resonances that will decay by annihila- 
tion into only Standard Model particles with strong or electromagnetic interactions. A particularly 
interesting example is the possibility of stoponium, a 2S+1 Lj = l So bound state of the lighter top 
squark and top antisquark, which we will denote rjf . If the light top squark (or stop) t\ has no kine- 
matically allowed flavor-preserving two-body decays, then it will live long enough to form hadronic 
bound states including stoponium. Specifically, this will occur if the decays t% — > bC\ and t\ — ► tN\ 
are kinematically forbidden, where C\ and N\ are the lighter chargino and the lightest neutralino 
respectively. As pointed out originally by Drees and Nojiri [l|, and re-examined recently in 
0, 0] , decays of s-wave stoponium to 77 can provide a viable signal at the Large Hadron Collider 
(LHC) with sufficient beam energy and integrated luminosity. The signal will appear as a peak 
with width determined by detector resolution on top of a smoothly falling background. This could 
provide a uniquely precise measurement of the top squark mass, which would in turn serve as a 



12] for other 



"standard candle" for kinematics of the supersymmetric sector. (See also refs. 
work related to stoponium at colliders.) 

As reviewed in more detail in refs. 0,0], there are at least two good motivations for considering 
top squarks light enough to form bound states. First, in models of "compressed supersymmetry" 



13], the top squark automatically comes out relatively light (compared to so-called mSUGRA 
models). These models can provide the right amount of dark matter provided that the mass 
difference between t% and N\ is less than the top quark mass, with m; between about 200 GeV and 
400 GeV. Second, minimal supersymmetric models that can provide baryogenesis at the electroweak 
scale [li], [3] require a quasi-stable stop with a mass that is currently estimated to be in the range 
from approximately 118 GeV to 135 GeV. 

The production and decay of stoponium as studied in refs. [l|, d, 0] at leading order (LO) 
are subject to large QCD radiative corrections and a strong dependence on the renormalization 
scale. The purpose of the present paper is to compute the corrections to stoponium production in 
hadronic collisions at next-to-leading order (NLO) in QCD, in order to help interpret limits on, or 
discovery of, a stoponium resonance at the LHC. This result will complement (and will make use 
of) our previous calculation of the QCD NLO corrections to stoponium decay [4|. 

In this paper, we will model the inclusive production of stoponium at a hadron collider as 
factored into the production of the free squark-antisquark pair with the same momentum in a 
color-singlet state and the binding of the pair into the stoponium bound state. For each pair 
of initial-state partons, the production amplitude for the squark-antisquark state with quantum 
numbers identical to the scalar ground state of stoponium is calculated in perturbative QCD, then 
related to the parton-level stoponium production cross section by integrating over the phase space 
that models the non-perturbative process of bound state formation using the wavefunction at the 



origin. In carrying out this program, we will closely follow the methods of ref. la ], which treated 



the analogous case of 1 Sq toponium (rj t ) production. (Ref. 16[] was written before the mass of the 
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top quark was known to be too large to allow it to form bound states, so toponium was a viable 
possibility at that time.) 

It is well known that the "static color singlet model" just described has failed spectacularly 
in explaining the observed J/i/j and ip(2S) prompt production in hadronic collisions. The color 
singlet model prediction for the prompt production cross section of the ip(2S) charmonium state is 
too small by a factor of about 50, while the prediction for prompt J/tp is too small by a factor of 
about 6, compared to the Tevatron data [17]]. (The T production cross section data 18|, 19|] seems 
to be acceptably fit by the color singlet model after inclusion of NNLO contributions [20j].) The 
failure of the static color singlet model in the charmonium case is due in part to the assumption 
that the hadronization of QQ can be factored into perturbative calculations of the production 
of open QQ with color and spin identical to the bound state and the nonrelativistic probability 
for annihilation at the origin. The discrepancy has provoked a variety of efforts to go beyond 



the static color singlet model in calculations; for reviews see [21[- 25]. For example, more general 
analyses performed using the NRQCD effective theory for nonrelativistic heavy quarks show that 
large enhancements to quarkonium production can come from the production of the constituent 
quarks in several different color and angular momentum states that can transition to the desired 
bound state. In this approach, the parton-level differential cross section for the production of a 
bound state $ in a collider is given by: 

da{ab^<5> + X) = ^ da{ab -► QQ[n] + X)(0*[n]), (1.1) 

n 

Here the differential cross section is factored into the short-distance perturbatively calculated 
cross sections initiated by partons a, b for the production of a quark-antiquark pair QQ in the 
color and spin state n, and the long-distance transition probabilities (O [n]) that describe the 
nonperturbative component of the hadronization. The sum over states n can be thought of as a 
perturbative series in the relative velocity v, which separates the long- and short-distance scales. 
So-called power counting rules are used to determine the relative order of each transition based 



on the effective theory interactions it represents [26|, |27(. The series converges more quickly as 
the mass of the quarks increases; for charmonium v is approximately 0.3 and for bottomonium 
v ~ 0.1. Transitions suppressed by the relative velocity can nevertheless have large cross sections 
and be much larger than the direct production assumed in the color singlet model, which may help 
explain the large observed prompt production cross sections. However, measurements of the J/ijj 



polarization 28f] at the Tevatron and other measurements 29( by PHENIX at RHIC do not seem 
to fit well the predictions of the NRQCD approach. 

Fortunately, however, the importance of corrections beyond the static color-singlet model should 
be relatively much less important for stoponium production, for several reasons. First, because the 
top squarks are much heavier, expansions in as should converge more quickly. Second, because 
r)i formed from scalars is a J = state, the LO cross section for the color singlet state does not 
require an extra gluon, as is the case for J/tp, ip{2S) or T production at leading order. Finally, in 
the stoponium case, all color octet states that can transition to the color singlet stoponium state 
1 S'o are of relative order v A . They are already negligible in the case of % production at next-to- 
leading order [30j], and since v scales with as, these color-octet transitions should be suppressed 
by an additional large factor for stoponium. These considerations justify our use of the color 
singlet model to calculate the radiative corrections to scalar stoponium production and decay in 
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perturbative QCD. Note that this essentially amounts to letting the wavefunction(s) at the origin 
R(0) in a potential model play the role of the spin-0 color-singlet matrix element (s) in eq. (jl.ip . 
while neglecting the effects of the higher spin or color matrix elements. A more accurate treatment 
including those effects may eventually be necessary, but is beyond the scope of the present paper 
and in any case would require some way of estimating the other relevant matrix elements, which 
are presently unavailable. 

The organization of the remainder of this paper is as follows. In section II we provide the NLO 
parton-level cross sections for stoponium production in hadronic collisions in the static color singlet 
model. In section III we discuss how to turn these calculations into hadron-collision cross sections. 
In section IV, we provide numerical results for the NLO stoponium production cross section in 
proton-proton collisions, studying the dependence on beam energies relevant for the CERN Large 
Hadron Collider (LHC), on the stoponium mass, and on the renormalization and factorization 
scale. We also review the next-to-leading order hadronic and diphoton branching ratios for the 
annihilation decays found in [4], which we use to estimate the NLO QCD cross section times 
branching ratio for the observable signal pp — ► rj^ + X — > 77 + X in scenarios where the hadronic 
partial width dominates the full width, and discuss corrections that apply to this idealized result 
in compressed supersymmetry and in models with electroweak-scale baryogenesis. 



II. PARTON-LEVEL CROSS SECTIONS 

In this section, we calculate the individual parton-level cross sections that contribute to the 
hadronic production of stoponium in NLO QCD. We will closely follow the method used by Kiihn 



and Mirkes to calculate radiative corrections to toponium production in ref. [la]. Although the 
virtual corrections to the leading-order diagrams for stoponium differ from the toponium result, we 
will show that the other radiative corrections contributing to the next-to-leading order cross section 
have the same form for stoponium as for their toponium counterparts, when both are written in 
terms of their corresponding leading order results. The relevant Feynman rules for our calculation 
can be found in Q]. Both ultraviolet and infrared divergences will be dealt with by dimensional 
regularization. 

In order to compute the parton-level cross-sections for stoponium production from the rate for 
open squark-antisquark production, integration over the final-state phase space must be restricted 
to the subspace containing the bound state. Let rfPjv(&l> k2',Pi,P2, ■ ■ ■ ,Pn) be the Lorentz-invariant 
differential phase space element for an V-body final state. In our case, N = 2 or 3, and k\ 
and k% will be the initial (massless) parton momenta and p\ and P2 are the final squark and 
antisquark momenta and ps is a possible final state (massless) parton momentum. To project onto 
the bound state phase space, we require the final-state squarks to have identical 4-momenta, letting 
p = pi = P2, so that 2p is the stoponium momentum, and use the radial wavefunction at the origin 
i?(0) (normalized so that i/j(r) = R(r) /\^4tt for an 5-wave state) to characterize the long-range 
behavior of the squark hadronization. As a consequence of this projection, the relevant 2-body 
and 3-body differential phase space factors in d = 4 — 2e dimensions are replaced by 

dV 2 - ^^-M 2 ), (2.1) 
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FIG. 1: The leading-order diagrams for the production of stoponium. 

where the parton-level Mandelstam variables are defined by 

a = (h + k 2 ) 2 , t = (2p-k 1 ) 2 , u=(2p-k 2 ) 2 , (2.3) 

and satisfy 

§ + i + u = 4p 2 = M 2 (2.4) 

where M = 2mj — -^binding ~ ls the bound state mass. Here and in the remainder of the paper, 
we will often write t without the subscript 1 and M instead of M v ~, for simplicity. 

A. Gluon fusion at leading order 

The parton-level diagrams for the leading order production of stoponium are in Fig. Q3 In 
this calculation, the squark and antisquark are restricted to color combinations appropriate for 
the singlet bound state using a wavefunction projection factor Sj/^/Nc, where i,j are SU(N C ) 
fundamental representation indices (with N c = 3 in the real world). The squared amplitude 
obtained from these diagrams, averaged over initial color and polarization and summed over final 
color, is 

1 V|M|2 87T 2 5|^ 

(d - 2) 2 (AT? - 1)2 \ J l l iv c (iv2 - 1)(1 - e ) ' y ' 

where [i is the regularization scale and as is the bare coupling. 

The leading order differential cross section for stop-antistop production in a color singlet state 
is therefore 

da LO (gg - tt*) = j- Hd _ 2)2 1 {N 2 _ 1)2 E l-^l 2 ^, (2.6) 



where s is the parton-level center-of-momentum energy squared. The leading-order cross-section 
for scalar stoponium production is obtained from this by the projection that restricts the squark 
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and antisquark to identical 4-momenta and includes the probability of annihilation at the origin, 
using eq. f|2. 1[) : 



&Lo(gg -> Vt) 



47r 2 a|/x 4e \R{0)\ 2 5(1 - z) 
N C (N? - 1)(1 - e) M 3 s 



(2.7) 



where 



(2.8) 



Note that the result (|2.7j) for LO stoponium production is a factor of 2(1 — 2e) smaller than the 
corresponding result for toponium obtained in eq. (10) of ref. 161 ]. 



B. Gluon fusion at NLO 

The corrections to leading-order gluon fusion can be divided into two parts - the virtual cor- 
rections coming from gluon loops in the leading-order diagrams and the corrections that involve 
the real emission of an additional gluon in the final state. The virtual corrections to stoponium 
production through gluon fusion are identical to the virtual corrections to stoponium annihilation 
to two gluons, which we have already calculated [4|. Adding up the two-particle cuts in Tables I 
and II of ref. [4j] , and combining with eq. (|2.7f) , we find 



OLO+virtual^ -> Vt) = a oS(l ~ z) 



2 \M 2 



N c b ( , 5vr' 



-IK. 



+ 1 + 

2em V 12 



Nr. 



3 + 



7T 



c F 



(2.9) 



where C F = (iV c 2 - 1)/2N C , and 



11 AT 2 

bo = yiV c - -n f , 



(2.10) 



with rif the number of quark flavors, and we define for future convenience 

4vr 2 a| \R(0)\ 2 



cr = 



N C (N 2 - 1)1 M 3 



(2.11) 



In eqs. (|2.9p and (|2.1ip . as(Q) is the MS coupling renormalized at the scale Q, related at one-loop 
order to the bare coupling as and the regularization scale /i by 



a s = a s 



1 _ as b ° 



7T 4 \ e uv 



1 



+ ln(4V/<2 2 ) " IE 



(2.12) 
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(a) (b) (c) 

FIG. 2: Non-zero diagrams related by crossing symmetry to the real gluon emission corrections to scalar 
stoponium production. Diagrams related to (a) and (b) by arrow reversal are also included with them but 
not shown. 



Note that we do not include the two-cut diagrams relating to the insertion of squark loops in the 
gluon propagator (diagrams ql-q4 of ref. 0]) since we will use from the rif = 5-quark effective 
theory to be consistent with the parton distribution functions we will use. We have also made the 
replacement In (m 2 /M 2 ) — > 1/ctr in the mass singularity arising from quark loops and the 1/v 
singularity is absorbed into the definition of the bound state wavefunction. 

The squared matrix element for real gluon emission can be obtained from that of the diagrams 
for ggg — > tt* in Figure [2] by crossing one gluon to the final state in all possible ways. The three 
external gluons carry labels i,j,n = 1,2,3 and carry momenta and polarizations fcj and etc. 
The stop and antistop each with mass m have the same momentum p (since we let the relative 
velocity go to zero) . The matrix element resulting from Figure [2] is 



+M {b) +M i - c) 



(2.13) 



where 



• ka (a) _ (fcj + fc n )^(A: n kj) v (kj + kj) p 
» vp ~ 2[(p - h) 2 - m 2 ][(p - k n ) 2 - m 2 ]' 



(2.14) 



. • / <(6) _ {9nvy^i kj) ' kn + (ki + 2kj) p{k n ) v (2/Cj + kj) l/ (k n )p]{ki + kj)p 

ilvp ~ [(p-k n ) 2 -m 2 }(k l + k 3 ) 2 ' [ ■ b) 

iM$ p = [g^{ki - kj) p + g up (ki + 2kj)^ - g P p{2k t + kj) u ]/(ki + %) 2 (2.16) 

correspond to the contributions from diagrams of types EJa), (b), and (c) respectively. The total 
matrix element that includes all diagrams is 

M=J2 M (i,J,n), (2-17) 

where Yl(i,j,n) represents a sum over the six permutations (i,j,n) = (1,2,3), (1,3,2), (2,1,3), 
(2,3,1), (3,1,2), (3,2,1). Taking the labels 1,2 to refer to the initial-state gluons and 3 for the 
final-state gluon, the spin-summed squared matrix element for gg — > tt*g is equal to Yl \-M\ 2 by 
crossing, when expressed in terms of the parton-level Mandelstam variables defined in eq. (|2.3p . 
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The sums over the physical polarizations of the gluons are performed using 

£ e,(X, k)e* u (X, k) = - 9lv + ^; + ;^ - ^ (2.18) 

A 

where k is the momentum of the gluon in question and n is an arbitrary 4-vector. This avoids the 
added complication of needing to cancel unphysical polarizations using additional diagrams that 
include ghost loops for each triple gluon vertex. Note that it is convenient for any given gluon 
polarization sum to choose any other of the massless gluon momenta for n. Including factors for 
the averaging over initial gluon polarization and color, this parton-level differential cross section is 

da(gg - m ) = Y H d-2) 2 (N 2 -l) 2 E l^l 2 (2.19) 

where dVj, is understood to be replaced according to eq. (|2.2p to obtain the stoponium + X 
differential phase space. The result is 



st + tu + us 



(s — M 2 )(t — M 2 ){u — M 2 ) 



M 8 + g 4 + g + u 4 
stu 



dt, 



(2.20) 



where we have dropped a term proportional to e that does not have a potentially singular denom- 
inator after angular integration in the e — > limit. 

Angular integration is carried out by replacing the Mandelstam variables with the dimensionless 
variables y and z defined on the interval [0, 1] by eq. (|2.8p and 



y = (l + cos0)/2, 



(2-21) 



where 6 is the angle in the center-of-momentum frame between the initial-state parton with mo- 
mentum k\ and the stoponium momentum direction. This implies that 



i=-M\l-y)(l-z)/z, 



u = -M 2 y(l - z)/z. 



(2.22) 



In terms of y and z, the partonic cross section is 



Hgg vw) = a o — N < 



vr \ M 2 J T(2-e) 

I - y(l - y)(l - z) 
[l-y(l-z)][z + y(l-z)] 



dyy e {l-y) 



-1-e 



+ l + (l-z) 4 [(l-y) 4 + y 4 ]), (2.23) 



where we have used the symmetry of the rest of the integrand under y — ► 1 — y to replace y 1 e (l — 
y)- 1 ^ with 2y- t (l - y)~ l ~ e . 
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To compute this integral, one uses plus distributions to simplify the integrand and isolate the 
soft and collinear divergences before integration is carried out. The plus distribution (F(x)) + of a 
function F(x) is defined by |32| 



(F(x)), = lim 



1-/3 

F(x)S(l - x - (3) - 5(1 - x - (5) I F{y)dy 

o 



(2.24) 



from which follows the identity 



G(x) (F(x)) + dx -- 

JO 



G(x) -G(l) F(x)dx 



(2.25) 



and the expansion identities 
y-(l- y)- 1 - = 



>-») + (rb 



Iny _ /^ln(l - y) 



✓ (1 -.)"'- ^(l-„ + ( T i- 



+ e 



i-y 

In z 
1 - z ' 



2e 



1-2/ 
ln(l - z 
1 - z 



(2.26) 
(2.27) 



Using these expansions in eq. (|2.23p and integrating over y using eq. (|2,25p . one obtains: 



°(gg Vt9) 



as, 



e z 



47T^ 2 \ /T e r(l + e) I 1 



7T" 



1 - e 



1 - z 



—8(l-z)--5(l-z) 
llz 6 + 2z 4 + 24z 3 + 23z 2 + 12 



2z 7 + 3z 6 + z 4 + 2z 3 + 5z 2 - 1 ■ 

+ z(l-z)(l + z) 3 Z 

[i - ,(i - r /Mi^) 

z(l-z) L V 7 V 



12z(l +z) 2 



1 



1 - z 
— lnz 



(2.28) 



Further simplification can be carried out by noting that if G(l) = 0, then eq. (12. 25ft implies 



G(z) (F(z)) + dz= f G(z)F(z) dx, 
o Jo 



(2.29) 



and therefore one can simply replace the plus distributions multiplying such functions accordingly, 
under the assumption that we will eventually integrate to an upper limit z = 1. The result of this 
simplification is 



aS AT 
IT 



4V\ V e r(i + e) I 1 



+ 



M 2 



1 - z 



1 - e 



5(1 



7T 



5(1 



+ Z(1-Z)~2 +Fgg(z) 



(2.30) 
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where 



llz 5 + llz 4 + 13z 3 



6z(l + z 



I9z 2 + 6z-12 /l 

r= + 4 - + z(l 



2 ) ln(l 



+4 



ln(l - z) 
1 -z 



+ 



2(z 3 - 2z 2 - 3z - 2)(z 3 - z + 2)zlnz 
(l + ,z) 3 (l-z) 2 



1 - z 



(2.31) 



This is identical in form to the corresponding real emission correction to the leading order parton- 
level toponium cross section [16], when both are written in terms of their respective LO results 

<T . 

To get the full next-to-leading order QCD cross section for gluon fusion, one must add the 
leading order plus virtual corrections from eq. f|2.9[) and real emission corrections from (|2.30p . 
resulting in 



Ac 



+(N C - 3C F ) (l + ^ | + N C F< 



gg( z ) ~ — p gg( z ) 



(2.32) 



where 



Pgg(z) = 2N C 



1 

- + 

z 



1 - z 



+ z(l-z) 



(2.33) 



Factorization requires the subtraction of 



a s /x 4 T(l + e) /4vr M 
^ l-e [ Of 



2\ e 



eiR 



Pgg(z) + C gg (z] 



(2.34) 



from eq. (|2.32p . where C gg {z) is scheme-dependent. In the MS factorization scheme that we will 
use for numerical work below, C gg {z) = 0, but it is a non-trivial function of z in other schemes 
such as DIS. Taking e — > gives the final parton-level cross section 



TT 



5(1 - z) 



P >rM{% i )+N c F gg (z)-C gg (z) 



\M 2 



(2.35) 



The form of this result is very similar to the corresponding one for toponium production, eq. (38) 
in ref. 16]], when both are written in terms of their respective LO contribution cro, differing only 
in the coefficient of aoCFCes/n- 
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FIG. 3: The non-zero diagrams corresponding to the leading-order production of stop-antistop in a color- 
singlet state through quark-gluon scattering. The additional diagram obtained from (a) by arrow reversal 
on the squark line is included but not shown here. 



C. Quark-gluon scattering 



Quark-gluon scattering contributes at relative order as compared to the LO gluon fusion cross 
section. The diagrams for the process qg — * qtt* are given in Fig. [3j Let k\ be the initial-state 
gluon momentum and ki (Aft) be the initial (final) momentum of the approximately massless quark, 
and let the squark and anti-squark each have momentum p. Then in terms of the parton-level 
Mandelstam variables s, t, and u defined in eq. (|2.3|) . the squared amplitude after summing over 
spin, polarization, and color is 



t 



2su 



\t-M 2 f 



1 + e 



(2.36) 



Inserting the required factors to average over the initial spin, polarization, and color, the parton- 
level differential cross section for qg — * rj^q is obtained by making the replacement of eq. (|2.2p for 
dVs in 



da(qg -> t#j) 



2s 2(d - 2)iV c (iV c 2 - 1) 



-£LX| 2 dP 3 . 



(2.37) 



The result is 



da(qg -> rftq) = cr ^C F 



k 



2vr r(2 - e) 



47T/U 2 s\ M 



lit 



St 



2su 



(s + uY 



1 + e 



dt. 



(2.38) 



The angular integration required to obtain the parton cross section is performed in exactly the same 
way for gluon fusion. Replacing the Mandelstam variables s, t, u with y and z using eqs. (|2.8j) and 
(pT22j) results in 



a s n ( 



l 

dyy~ e (l -y) 



•k 



-z 1+e (l-z) 
r(2-e) V ' 

_ x _ e f l + y\l- Z f 
[l-y{l- Z )f 



-2c 



(2.39) 
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FIG. 4: Non-zero diagrams corresponding to the leading-order production of stop-antistop in a color-singlet 
state through quark-antiquark annihilation. Note that these are related by crossing to the diagrams of figure 
[3l The additional diagram obtained from (a) by arrow reversal on the squark line is included but not shown 
here. 



Using the expansions in equations (|2.26p and (|2.27p . we integrate this to obtain 



v x as /VVA(l + e) 



l + (l-z) 2 ( 1 



+ 21n(l-z) 



+2 + z - 2/z - zlnz 



The infrared divergence is removed in factorizing the cross section by subtracting 



(To; 



as/i 4 T(l + e) /4vr/x : 



,2\ e 



Ql 



-P gq (z) + C gq (z) 



2tt 1 - e 

from eq. (|2.40|) . where the splitting function P gq {z) is defined by 

P gq (z) =C F [l + (l-z) 2 }/z 



(2.40) 



(2.41) 



(2.42) 



and Cg q (z) = in the MS factorization scheme that we will use for numerical work. Now taking 
e — ► 0, we arrive at the final result for the order-a| quark-gluon parton cross section 

HQ9 -» m) = ^ [ p m( z ) ln ( m2 (! " z ?/Q 2 f) + Cf(2 + z - 2/z - zlnz) - C gq {z)] . (2.43) 

Just as we found for the real gluon emission corrections to the leading order diagrams, this result is 
identical to the corresponding result for the production of quarkonium (eq. (47) of ref. lGf] ) , when 
both are written in terms of their respective LO results <to for gluon fusion. The cross section for 
antiquark-gluon scattering (qg — ► -q^q) is also given by eq. (|2.43[) . 



D. Quark-antiquark annihilation 

The production of color singlet stoponium through quark-antiquark annihilation is impossible 
at leading order. However, the process can occur through the emission of a real gluon, through the 
diagrams of figure [H which are related to the diagrams for quark-gluon scattering in figure [3] by 
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crossing the initial gluon line and the final quark line. Therefore, the quark-antiquark annihilation 
squared amplitude Yl l-M] 2 can be obtained directly from eq. (|2.36p by the substitutions (s, t, u) — > 
(u, s,t), with an extra factor of —1 because we are crossing a fermion. This immediately yields 



Ag e C F 



2tu 



1 + e 



(s-M 2 ) 2 

There is no possible infrared divergence, so we can immediately take e = 0. Using 

1 1 



(2.44) 



da(qq -> r^g) 



Y,\M\ 2 dP 3 



2s 4iV 2 

and applying the substitutions in eqs. (|2.2p and (|2.22p . we have 



(2.45) 



- »7#) = ^o— C 2 F z(l -z) [1 - 2y(l - y)} dy = a ^C 2 F -z(l - z). (2.46) 

7T Jo ' 7T 3 

This has the same form as the corresponding parton cross-section for toponium (eq. (49) of ref. la]), 
when both are written in terms of their respective leading order gluon fusion results co- 



lli. CROSS SECTIONS AT HADRON COLLIDERS 

The parton model and QCD factorization relate the hadron-level^ cross section to the parton 
cross sections through: 




where fa(x a ,Q F ) is the parton distribution function for parton a carrying momentum fraction x a 
in a proton probed at a scale Qp. Since the parton cross sections are integrated in terms of z, we 
make the change of variables x a = x, x& = t/xz, where r = M 2 /s, where y/s is the hadron-level 
center-of-momentum energy. This gives 

<?(PP -^Vi + X) = J2 [ dz I dx ^2 Q F )f!(r/xz, Ql) a(ab m + X). (3.2) 

a,b ' 

[It is sometimes useful in numerical work to decouple the limits of integration by making the further 
change of variables x = 1 — v(l — t/z).] 

Equation (|3.2|) is applied using eqs. (|2.35p . (|2.43[> for both quark-gluon and antiquark-gluon 
scattering, and (|2.46p . Since the lower limit of z integration is r, the plus distributions in the 
parton-level cross sections must be shifted in order to be easily integrated. Define the "r+" 



We will write formulas as they would apply to pp colliders such as the LHC, but simple changes can be applied in 
the obvious way to obtain pp collider results. 
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distribution as 



lim 

13^0 



1-/3 

F(x)0(l-x- P)-5(l-x- (3) I F(y)dy 



(3.3) 



which is related to the plus distribution by 



(F(x)) T+ = (F(x)) + + lim 6(1 -x- 13) F(y) dy. 



(3.4) 



Now the plus distributions occurring in the gluon- fusion cross sections [see eqs. (|2.3ip . (|2.33|) and 
(|2.35p ] are replaced using 



1 



1 - z 
ln(l - z) 
1-z 



1 



1 - z 



+ ln(l - t)6(1 - z), 



T+ 

ln(l - z) 



1 



, +l]n 2 (l-r)6(l-z), 

1 Z I T+ L 



(3.5) 
(3.6) 



and integration of terms containing r+ distributions can proceed using 



(F(z)) T+ G(z) dz = / G(z) - G(l) F{z) dz 



(3.7) 



Note that the simplification in equation (|2.29|) is not affected by changing the lower limit of inte- 
gration over the plus distributions. 



IV. NUMERICAL RESULTS 

In this section we present our numerical results for the stoponium cross section at next-to-leading 
order in QCD in pp collisions at energies relevant for the LHC. (We note in passing that the QCD 
cross section times branching ratio for pp — > r] j + X — > 77 + X at y/s = 1.96 TeV is much less than 
1 fb even for stoponium as light as 200 GeV, so searches for rjf — > 77 at the Fermilab Tevatron are 
presumably hopeless unless there is some very strong new non-QCD production mechanism.) 



A. Stoponium wavefunction effects 

In the static color singlet model all of the nonperturbative information about the formation of 
the bound state is contained in the amplitude of the wavefunction at the origin (and its derivatives, 
for non-zero angular momentum states) . In section [III the parton-level cross sections were obtained 
in terms of |i?(0)| 2 /M 3 , a quantity that can be calculated approximately by using potential models 
that simulate the effects of QCD. A naive Coulombic model for the QCD binding force would imply 
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that for the S-wave level n state (n = 1 is the ground state), 

(4.1) 



\R nS {0)\ 2 _ 4a| 



M(ns) 3 27n 3 ' 

but asymptotic freedom and other perturbative and non-perturbative QCD effects make this ap- 
proximation quite crude. We will instead obtain the value of the wavefunction at the origin by 



using a parameterization 3l( based on a potential extrapolated from the study of charmonium 
and bottomonium. The numerical results for \R n s(0)\ 2 /M(ns) 3 in ref. 3lJ depend on the 4-flavor 



QCD scale A§£ D as an input. For the range Og\mz) = 0.1185 db 0.0025, the relevant range is 
Aqqp_~ 300 MeV ± 40 MeV, and so we will use the Aq^, d = 300 MeV parameterizations from 



ref. 3l| in our numerical analysis. 

While the largest stoponium production is for the 15 state, it can also be produced in excited 
states, and the signal will be enhanced by their decays either directly to two photons or to lower 
states which then decay to two photons [3]. As mentioned in this reference, the production of 
higher angular momentum states is essentially of relative order a| 3l|], so we do not consider 
their direct production as a correction. However the contributions of higher S-wave states will be 
non-negligible. Although their annihilation decay branching ratios will be the same as the ground 
state, they may also decay to P-wave bound states that cannot annihilate directly to a two-photon 
final state, but can decay to the ground state or other S'-wave states. The total diphoton signal 
is therefore presumably bounded above by the sum over the production cross sections for all nS 
states, but it is unknown how much the higher n states will contribute to the diphoton signal. 
However, since the phase space for 25 — > IP decays is highly suppressed due to a very small mass 
difference, the diphoton decay branching ratios of the 15 and 25 states should be nearly identical. 
In our numerical results we will therefore conservatively assume that the relevant production is 
due to these two states, and use in the results of section |H] the effective wavefunction at the origin 
factor: 

TO)! 2 y l^s(0)| 2 U o) 
M 3 ^ MinSf K ' ' 

n=l v ' 

neglecting the contributions from n > 3. The higher energy states may also contribute significantly, 
but with diphoton branching ratios suppressed by an unknown factor due to available decays to 
P-wave bound states that do not eventually decay to diphotons. Furthermore, it is quite possible 
that the potential models become less reliable for the higher excited states. 



In Figure [5] we plot the effective wavefunction at the origin factor given in 31[ as a function 
of the IS bound state mass. Three cases are shown in this plot: the contribution from the IS 
bound state alone, the sum of the IS and 2S bound states that we will use in numerical work, 
and the sum of the first 10 5-wave states. In each case, the 40 MeV variation in the QCD scale 
Aqq D (obtained by interpolation from the values given in ref. 31]) is also plotted on the graph. 
Although the probability for producing higher-energy bound states drops sharply with higher n, 
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FIG. 5: The bound state wavefunction at the origin 
factor \R(0)\ 2 /M 3 , as a function of the IS bound 
state mass, as parameterized in [311 ] . The lower 
three lines are for the IS wavefunction only, the 
middle group is the sum of the IS and 2S wavefunc- 
tions, and the top group is the sum of the lowest 
10 iS-wave bound state wavefunctions. The solid 
lines are the Aq^ d = 300 MeV parameterization 
and the dashed lines show the ±40 MeV variations 
in this value. 
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FIG. 6: Binding energies 2toj — M^^nS) of the nS 
bound states of mass M Vi (nS) as a function of the 
constituent squark mass nif. The top line is the 
IS state, the second is the 2S state, and so on to 
the 10S state at the bottom. This plot uses the 
Aq^ d = 300 MeV parameterization in ref . 31 1 . 
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their sum can add up to a significant enhancement to the overall production rate.^ It is not unlikely 
that we have underestimated the cross section for stoponium production and the rate of diphoton 
annihilation decays, but precise calculation of the additional enhancements will have to await a 
better understanding of stoponium spectroscopy. 

In Figure [6] we plot the binding energies of the first 10 5-wave bound states, computed using 
the parameterization in 31]. These are about 1% of the IS bound state mass and significantly 
less for any other bound state. These differences are within the scale dependence of our result, 



It is interesting to note that the results for the wavefunction at the origin obtained in [3lJ are considerably 
smaller than would be obtained from the naive Coulombic formula (|4.1[) using as(Q) evaluated at Q — l/{r„s) = 
2as(Q)M rl{ /9n 2 . However, the results of eq. (|4.1|l depend very sensitively on this somewhat arbitrary choice of 
scale. Also, the contribution from higher n states from [3"fl ] is larger relative to the ground state contributions 
than is suggested by the Coulombic formula. 
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FIG. 7: The scale dependence of the LO (dashed) and NLO (solid) pp — ► m 
GeV (left) and 450 GeV (right), with the factorization scale Qf set equal to the renormalization scale Q 
for V« = 14 and 10 TeV. 



and more importantly they will be below the energy resolution of the detectors. Therefore we will 
ignore the small binding energies in the calculation of the parton cross sections, taking M = 2m^ 
everywhere and ignoring any small differences in excited state masses. 



B. Numerical results for stoponium production 

To integrate numerically the parton-level cross sections in equations (|2.35p . (|2.43p . and (12.461) 
using eq. ([33]) . we will use the Martin-Stirling-Thorne-Watt (MSTW) 2008 NLO parton distribu- 
tion functions 33J, using their global best fits. These PDFs use the MS subtraction scheme, and 



so we set C gg (z) = in eq. (p35j) and C gq {z) = in eq. (pQ3j) . To be consistent with MSTW 



3J|, we run as(Q) with two-loop beta function in the nf = 5 effective theory, starting from 
a s {mz) = 0.12018. 

In Fig. [7] we plot the renormalization scale dependence of the NLO pp — > rjf + X cross section 
for two different masses M m = 240 and 450 GeV, varying 50 GeV < Q < 480 GeV and 100 GeV < 
Q < 900 GeV respectively as reasonable ranges for the choice of the common renormalization and 
factorization scale Q = Qf- The figure shows the results for both yfs = 10 TeV and 14 TeV proton- 
proton colliders. These results show a significantly improved renormalization scale dependence for 
the NLO result compared to the LO cross section in each case. 

We present the individual parton-induced contributions relative to the full hadronic cross section 
in Fig. [U for yfs = 14 TeV proton-proton colliders and for the same masses and over the same 
ranges of Q as in figure [7] Not surprisingly, for Q = Qf = M^., the relative contributions from 



each parton-level process are quite similar to those found for the toponium cross section in ref. 16f] . 
It is also interesting to note that for the scale choice Q = Qf ~ M V{ /3, the NLO gg, qg + qg and 
qq corrections are all simultaneously small. We have checked that this also holds for lower beam 
energies. 
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FIG. 8: Relative contributions from each parton-level process relative to the total NLO pp — ► rn + X for 
M Vi = 240 GeV (left) and 450 GeV (right), for y/s = 14 TeV, as a function of varying Q = Q F , as in figure 
[7J The LO contribution is the order a 2 s result, the gg contribution is due to real and virtual corrections to 
the leading order diagrams (the remaining part of eq. (|2.35p . the qg + qg contribution is from quark-gluon 
and antiquark-gluon scattering eq. (|2.43[) . and the qq contribution is from quark- antiquark annihilation 
eq. (Qg) . 
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FIG. 9: Results for NLO order cross sections for pp — > 
rii+X with M v - between 200 and 800 GeV. Here we have 
set the factorization and renormalization scales equal to 
the stoponium mass, Q = Qf — M m . 
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In Fig. [_U we show the NLO pp — » ?^ + X cross section as a function of , with the renor- 
malization and factorization scales set equal to the stoponium mass. At this writing, the ultimate 
performance of the LHC is the subject of considerable speculation, so we show results for four 
different beam energies. Figure [T(_| shows the corresponding K-factor, defined as the ratio of the 
next-to-leading order to the leading-order cross section, again with Q = Qf = M^-. Our results 
show that the enhancement in the stoponium production cross section due to the NLO corrections 
with this scale choice is between 30% and 50%, depending on the energy of LHC, with larger 
enhancements for larger M v . and for smaller y/s. 
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FIG. 10: The K-factor, defined as the ratio of the 
NLO to LO cross sections for pp — > rji + X, as a 
function of the stoponium mass M v -, computed with 
Q = Qf = M V{ . Results are shown for y/s = 14, 12, 
10, and 7 TeV. 
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C. Numerical results for the diphoton final state 

The observable signal for stoponium production is the diphoton peak produced by the decay 
Vt ~* 77- The cross section for this signal is the product of the production cross section we have 
calculated above and the model-dependent branching ratio for diphoton decays of stoponium: 



a(pp -> r}~ t -> 77) = a(pp -> r? f )BR(ry t - -> 77) 



(4.3) 



Unfortunately, the 77 branching ratio is in general highly dependent on the parameter space of 
soft supersymmetry breaking. However, one can proceed by first considering the idealized model- 
independent case in which the total width is dominated by the gluon-gluon partial width. Then 
we can approximate the full width by the NLO width to hadrons (gg + X or qq + X), and calculate 
an approximate branching ratio: 



BR(r? t ~ -» 77) 



rd)fa-^ 77 ) 

r( 1 )(77| — ► hadrons) 



(4.4) 



where — > hadrons) and r^(^ — » 77) are the NLO partial widths. In our earlier paper [4j, 
we found 



9cr| 



7T 



^ln(«i | + 
2 \M 2 



24 



IS ,^+( T -2-21n(2))C F 



+ I l( n f + nt) + n t h(im 2 t /M 2 ) + i ln(2) j 



(4.5) 



in which M is the bound state mass, nj = 5, nt is zero or one depending on whether or not the 
mass of stoponium is large enough for top-antitop decays to be kinematically allowed, and the 
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FIG. 11: The NLO branching ratio of stoponium 
to two photons in the idealized approximation 
that the hadronic partial width dominates the full 
width, from eq. (|4.5p . We have used as{mz) = 
0.12018 and set the renormalization scale equal to 
the stoponium mass. 



FIG. 12: The NLO cross section times branching 
ratio for pp — > r\i + X followed by r\i — > 77, in the 
idealized case that the stoponium decay width is 
dominated by the decays to gluons. The factoriza- 
tion and renormalization scales are set equal to the 
stoponium mass. 
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function h(r) is defined by 

h(r] 



9 V 



-)VT 



r 



_l n (l + yr^7) + -ln(2). 



(4.6) 



The resulting idealized diphoton branching ratio is shown in Fig. [TTJ (Here we have differed 
slightly from figure 8 of ref. [4(] by using values of as(Q) that follow from as(mz) = 0.12018 to 
agree with those used in MSTW 3j], using a slightly different two- loop running that includes the 
stop squark and includes the top quark only if the stop mass is greater than the top mass. We 
have also fixed the renormalization scale equal to the stoponium mass; the NLO scale variation is 
quite small as shown in U) 

Now combining the NLO production cross section of fig. [9] with the NLO branching ratio of 
fig. [TT] we find the idealized NLO pp — > % + X — » 77 + X cross section shown in fig. [12j In this 
paper, we will not attempt a new study of the diphoton signal viability including backgrounds and 
after cuts. However, it is interesting to compare the corrected results we have found in the present 
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paper and in |4J to the corresponding LO results that were used as the basis for the study in 
ref. [jj. There are several new effects taken into account in this paper. First, the NLO corrections 
to the diphoton branching ratio lead to a decrease of 30% to 35% when evaluated at Q = M^.. 
Second, and counteracting this, the NLO corrections to the production cross section lead to an 
increase as shown in fig. [TU1 Third, we are including here the effects of the 2S state as discussed in 
section HV Al Also, ref. [3|] used a different set of parton distribution functions and a different as- 
The net result of these effects is shown in Table HI which compares the cross section (before cuts) 
used in the study ref. to that found here. The results found here for the cross section times 



M Vf [GeV] 


a x BR (fb, NLO) 


a x BR (fb, LO in ref. [3]) 


ratio 


235 


22.4 


23.1 


0.969 


240 


20.5 


21.0 


0.974 


255 


15.8 


16.0 


0.989 


270 


12.3 


12.3 


1.00 


400 


2.14 


1.94 


1.10 


450 


1.24 


1.09 


1.13 


500 


0.750 


0.647 


1.16 


600 


0.307 


0.255 


1.20 


700 


0.141 


0.113 


1.24 


800 


0.0700 


0.0549 


1.27 



TABLE I: Comparison of cross section times branching ratio for pp — > r\~ t + X — » 77 + X in proton-proton 
collisions at s/s = 14 TcV, before cuts and in the idealized case that the hadronic partial width dominates 
the decay width of stoponium. The second column shows the results in fig. [TJl and the third column shows 
the LO results as used in the study of the signal including cuts and LO backgrounds in ref. [|| . 

branching ratio are about equal to those used in 0] for 235 GeV < M Vi < 270 GeV relevant for 
electroweak-scale baryogenesis models, and are somewhat larger for larger M v . > 400 relevant for 
compressed supersymmetry models. 

In the two motivated scenarios mentioned in the Introduction, there are somewhat different 
expectations for the reduction of the pp — ► X — ► 77 + V cross section compared to the idealized 
case just mentioned. For compressed supersymmetry [131 ], in which enhancements to dark matter 
annihilation are mediated by t-channel stop squark exchange, the ground state stoponium will 
have a mass between about 400 GeV and 800 GeV, tending towards the former in models with less 
fine-tuning. The stoponium total width is indeed dominated by the hadronic partial width, which 
has a branching ratio that is typically more than 80% and usually of order 90% or even higher Q]. 
The resulting diphoton cross sections can therefore be estimated from figure [T2l simply by applying 
this correction factor as appropriate. 

It is somewhat harder to make a prediction for the diphoton signal in models of electroweak- 
scale baryogenesis , [3], because in those models the h°h° decay width typically dominates if 
it is kinematically allowed, and also because the WW and ZZ modes are more significant. We 
have calculated the radiative corrections to h°h° decay channel in ref. 0], but the weak vector 
boson decays are also significant (contributing up to about 30% of the total width for some choices 
of parameters). But to make a rough general statement about the cross section for diphoton 
annihilations, we can assume that the corrections to the weak decay channels should probably not 
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be much larger than those for the hadronic decays; then the 30% reduction in over the leading 
order ratio leads us to conclude that a corresponding reduction in BR(?7j — > 77) is reasonable. 
These assumptions and the results in refs. 0, Q] allow us to make a conservative estimate for the 
next-to-leading order branching ratio in the electroweak baryogenesis models: 

f 1.5 x 10~ 3 (M < 2m h o), 
BRfe - 77) > , (4.7) 

[ 5.3 x 10- 4 (M > 2m h o), 



Note that in these models, one expects M^- between about 235 and 270 GeV, which happens to be 
in just the range where the decays to h° might become kinematically allowed. The branching ratio 
is highly dependent on the Higgs scalar boson mass because of its large partial width immediately 
above the threshold for Higgs pair production. Our resulting estimates for the stoponium diphoton 
annihilation cross section in the Electroweak Baryogenesis scenario are given in Table HH 



M v . [GeV] 


a{pp -> rji) [pb] 


<7 77 , M < 2m h o [fb] 


<7 77 , M > 2m h o [fb] 


235 


8.25 


12 


4.4 


255 


5.63 


8.4 


3.0 


270 


4.30 


6.5 


2.8 



TABLE II: Estimates of the NLO cross section for the production and subsequent diphoton annihilation 
decay of stoponium in a t/s = 14 TeV pp collider, as a function of the stoponium mass in the electroweak 
baryogenesis model scenario. We have used approximate 77 branching ratios from eqs. (|4.7[) . 



V. OUTLOOK 

In this paper, we have derived the NLO QCD corrections to the cross section for stoponium 
productions in hadronic collisions relevant for the LHC. These could be important in understanding 
an eventual discovery of, or limits on, stoponium through its diphoton annihilation decays. We 
found that when calculated at a scale Q = M V; , the corrections to the cross section are large and 
positive, partly canceling the large negative corrections in the diphoton branching ratio found in 
ref. [4]. We also found that the NLO corrections to the cross section are small when evaluated 
at Q ~ M^/3. For stoponium of mass 240 GeV (450 GeV), we found that reducing the center 
of mass energy from 14 TeV to 12 TeV will reduce the cross section by about 25% (30%), while 
reducing from 14 to 10 TeV will reduce the cross section by about 52% (55%), with of course 
larger reductions for heavier masses. We have not attempted here a detailed study improving on 
the estimate of the reach of the LHC based on LO signal and background in [jj. One prerequisite 
for such an improvement is a better understanding of the Standard Model and detector diphoton 
backgrounds. There has been considerable study of the Standard Model high energy diphoton 



production at the LHC 35J]- 38[], but focused on lower invariant masses relevant to a light Higgs 
scalar boson search. Fortunately, in the future the LHC will provide its own background estimate 
in the form of a sideband analysis, and this should be at least as robust as any calculation can be, 
at least for the purposes of bump-hunting. 

The largest uncertainty in estimates of the stoponium cross section is clearly our lack of detailed 
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understanding of stoponium energy levels and bound state matrix elements. In this paper, we used 



the parameterization of [3l[ for the wavefunction at the origin, but it must be recognized that this 
was obtained in part by a significant extrapolation from charmonium and bottomonium data. Even 
at the energy scales relevant for stoponium, the QCD binding potential is far from the semi-classical 
Coulomb form, and more work is needed to understand the relevant matrix elements. We have 
tried to remain relatively conservative by not including the possible effects from S'-wave resonance 
production at the 3S level and beyond. An eventual discovery of stoponium in the diphoton channel 
will not only provide an accurate determination of the top squark mass, but an opportunity to use 
data to learn more about a QCD bound-state system in an energy range where calculations are 
under relatively better control than the known charmonium and bottomonium states. 
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